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Abstract 

We present numerical simulations of a model of cellulose consisting of long stiff 
rods, representing cellulose microfibrils, connected by stretchable crosslinks, repre- 
senting xyloglucan molecules, hydrogen bonded to the microfibrils. Within a broad 
range of temperature the competing interactions in the resulting network give rise 
to a slow glassy dynamics. In particular, the structural relaxation described by ori- 
entational correlation functions shows a logarithmic time dependence. The glassy 
dynamics is found to be due to the frustration introduced by the network of xy- 
loglucan molecules. Weakening of interactions between rod and xyloglucan molecules 
results in a more marked reorientation of cellulose microfibrils, suggesting a possible 
mechanism to modify the dynamics of the plant cell wall. 
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1 Introduction 

Cellulose is the most abundant organic material in nature and it is the basic 
structural component of the plant cell walls. In the higher plants, cellulose crys- 
talline microfibrils, 5-15 nm wide and 1-5 /im long, are embedded in a matrix 
of hemicellulose, pectin polysaccharides and proteins [I]. Bacterial cellulose 
composite, produced by Acetobacter xylinum growing in a medium containing 
xyloglucan (hemicellulose), is a pure cellulose-xyloglucan network with degree 
of crystallinity around 75% [2115] . Xyloglucan molecules are bound by hydrogen 
bonds to the amorphous surface of different cellulose fibrils creating crosslinks 
between them, which imparts elastic properties to the network. 
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The dielectric properties of bacterial cellulose were studied in the frequency 
range 100 Hz-1 MHz and at temperatures between 100 and 440 K [4J. A broad 
peak of dielectric absorption was found for T < 260 K, indicating a wide dis- 
tribution of relaxation times. Creep experiments on cell walls of higher plants 
and on composite material from Acetobacter xylinum have shown a logarithmic 
time dependence of elongation [5]. Moreover, the authors have shown that, in 
presence of the cell wall loosening lipid transfer proteins (LTP), the elonga- 
tion remains logarithmic but with a higher creep rate. These slowing down 
of relaxation and logarithmic time evolution are dynamical features typical of 
glassy systems [B"fT6 8.9j. Evidence of glassy behaviour in biologically relevant 
materials may be of general importance. Schrodinger in his seminal book [10] 
described the building blocks of life as "aperiodic crystals" , orderly structures 
withdrawn from the disorder of heat motion. Much later, the general concept 
of glassy state has been developed [6llll|ll2j ) as a state with infinitely many 
local metastable configurations. Such a glassy state is prone to minor changes 
of environment since even small external actions can transfer the system as a 
whole into a different metastable state, while never reaching thermodynamic 
equilibrium. Interestingly, the transition to a glassy state found in saccharides 
at low temperature [Sj has been suggested to be the key factor in the pro- 
tection of biological tissue against freezing. For these reasons, it is plausible 
to assume that living organisms should make use of substances with a broad 
variety of similar structures depending on external conditions, as offered by 
glassy materials. Evidence of slow glassy dynamics of the cytoskeleton of liv- 
ing cells was recently observed experimentally [1511161117] . The simplest and, 
at the same time, one of the most important organic materials, cellulose, is a 
good starting point to investigate this possible connection theoretically. 

The physical mechanism which is responsible for glass formation is the pres- 
ence of frustration, i.e. the incompatibility between the locally preferred order 
and global constrains. In other words, the energy of the system cannot be 
minimized by optimizing all local atomic interactions [T2lfT5] . 

In the cellulose-xyloglucan system, frustration gives rise to defects in the mi- 
crofibril organisation due to long-range xyloglucan-xyloglucan interactions. 
To simulate this phenomenon, we propose a model that describes cellulose 
microfibrils and xyloglucan molecules at the mesoscopic, coarse grained, level. 
Correlation functions of the microfibril orientation calculated within this model 
display the slow logarithmic behaviour typical of glasses over two time decades. 
Similar logarithmic decay of correlators has been observed in Monte Carlo 
simulations for spin-glasses [6J, colloids and micellar particles [9"|18|19j . and 
overcooled liquids [20] ■ This slow dynamics is often explained in terms of 
mode-coupling theory 08]. 

In the following we present the model and the results of molecular dynamics 
(MD) simulations for the cellulose network, analyze the structure and dy- 
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namics identifying the glassy behaviour, and point out the crucial role of the 
xyloglucan network for the slow dynamics of cellulose microfibrils. 



2 Model and Simulation Technique 



We propose the coarse grained atomistic model illustrated in Fig. [T] to simu- 
late bacterial cellulose pQ, a prototype of cellulose in higher plants. The model 
consists of long stiff rods, each constituted by seven beads, representing cellu- 
lose microfibrils, connected by a network of stretchable crosslinks, representing 
xyloglucan molecules. 

In our coarse-grained model, there are four types of interactions between: 
1) beads within a rod (b-b), 2) beads in different rods (r-r), 3) xyloglucan 
particles and rod beads (x-r) and 4) xyloglucan particles (x-x). We have cho- 
sen these interactions with the following criteria: i) the x-r bonding sets the 
scale of energy to that typical of a hydrogen bond, ii) the rods have to be 
stiff, impenetrable, very long compared to typical bonding lenghts, and the 
interactions between them are dominated by the interactions via xyloglucan 
particles, iii) the xyloglucan has to form a disordered network connecting dif- 
ferent rods with the possibility of bond stretching, breaking and formation. 
This is less straightforward than the previous steps and will be described in 
some details later on. 

We describe most interactions by a Lennard- Jones potential of the form: 
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where s = {b-b, x-r, r-r, x-x} labels the different interactions. The minimum 
of the potential, -e s , occurs at r eq s = 2 1 / 6 cr s . The constant C s is chosen in such 
a way that the potential vanishes at r = i? CjS , where R CtS is a cut-off radius 
and £/lj, s — for r > R C)S . The parameters are given in Tabled] in units of the 
characteristic energy and length of x-r interactions, e = e x _ r and a = o~ x _ r . For 
bacterial cellulose, these units can be assumed as e ~ 20 kcal/mol ~ 10 4 K, 
which describes a strong hydrogen x-r bonding [21], and a ~ 15 nm that gives 
the typical radius of the microfibrils. In these units, the equilibrium length of 
each rod corresponds to the typical length of cellulose microfibrils, L = 1 //m 
with an equilibrium distance between them h = 67 nm (Fig. [1]). Moreover, the 
cut-off radius of £/Lj, x -r is chosen in such a way that each xyloglucan particle 
interacts with no more than one rod -R CjX _ r — | r eq,r-r (Table [1]). 

In order to keep the rods rigid, we add to t/Lj,b-b a simple bending potential 
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between the nearest beads within a rod given as 

11 A/1 i \ ( r «-l — r i) r t+l ~~ r i) / n \ 
Ub = A(l + COSa), COSOj = ; ; r, (2) 

where A ^> 1 (Tabled]), r$ is the position of the bead i, with i = 2-6. This 
strong potential allows to decrease the cut-off radius R c ,b-b of £/Lj,b-b an d lower 
the computational time. 

We come now to the interactions between xyloglucan particles, which play 
the role of crosslinks between rods. To obtain a connected network of rods 
with effective long-range interactions, we require each xyloglucan particle to 
interact with two others. These triplets of connected particles are chosen once 
and for all at the beginning of the simulation as described in the following. 

We place three x-particles per rod at the equilibrium distance of £/Lj, x -r from 
the second, fourth and sixth bead of each rod. For each x-particle we choose a 
first partner as the farthest x-particle within a radius L that is not yet assigned 
as partner of another particle. We repeat this procedure to assign the second 
partner. Once all triplets have been assigned, the distances between them in 
the initial structure are defined as equilibrium distance r eq for the additional 
interaction between particles within a triplet described by a Morse potential 

U M = D(l~e- b ^- r ^)\ (3) 

with parameters D, b specified in Table [U This means that every pair inter- 
acts through a potential with minimum at a different value of r eq . As shown 
in Fig. [21 the initial r eq -distribution contains seven sharp peaks. We further 
assume that a bond between x-x pairs can be broken if Um ^ 0.9D, which 
introduces a cut-off radius R^ 1 for x-x interactions 

dM In(l-VO^) 2.97 

The resulting network is very robust and less than 6% of the x-x bonds break 
at all studied temperatures. The long-range interactions between xyloglucan 
particles introduce frustration in the orientation of the rods, that in their 
absence, would energetically prefer to form a crystal of parallel rods. Although 
there is a large arbitrariness in our construction of the interactions, we can 
assume that frustration due to competing interactions is a quite general feature 
of non crystalline multicomponent systems, (like most living matter). 

Simulations of glasses have to be conducted possibly over several decades 
in time. However this fact, together with the complexity of the interatomic 
potentials, limits the size of the system. In our simulation we use iV ro d = 288 
rods, located in parallelepiped with 6 unit cells in the x and y directions and 
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14 in the z direction. In each unit cell, defined by vectors 

i= (6(7,0,0), j = (0,6(7,0), k= (0,0,11(7), 

we put four beads at (0.25, 0.25, 0.25), (0.75, 0.75, 0.25), (0.25, 0.75, 0.75), (0.75, 0.25, 0.75). 

We form the rods by letting the beads interact with each other in groups of 

seven along the z direction, which implies two layers of densely packed rods 

along the z direction. In Fig. [3^ we show a typical snapshot of the equilibrated 

system where the two layers of rods and the x-r bonds are visible, whereas in 

Figure [3b we show only the bonds between a few xyloglucan triplets. 

In essence, our system is like a liquid crystal with the unit vector of rod orien- 
tation n (Fig.[T]), playing the role of order parameter. The constructed configu- 
ration, composed of N = 2880 particles, was equilibrated at high temperature 
T = Q.lQs/kB in the microcanonical NVE ensemble. We integrated the equa- 
tions of motion rescaled by the mass of the rod bead taken to be m = 10~ 22 kg, 
which is twice the mass of xyloglucan particles, by means of the velocity Ver- 
let algorithm. The time step in units of r = \Jmo 2 /e ~ 0.5 ns was chosen 
At = 0.032, the shortest period of vibrations between x-r particles being 
30 At. Then, the system was quenched with rate 3- 10 9 K/s and studied in the 
NVT ensemble, and studied at T = (3.0,1.9,0.9,0.45,0.23,0.11) • lO- 2 e/k B , 
using the Nose-Hoover chain thermostat [22j. The masses of the thermostat 
chains were defined by Qi = NT/u 2 and Q 2 = T/oj 2 with to = 6.7, the char- 
acteristic frequency of x-r vibrations. For convenience, in the following, we 
give the values of temperature in reduced units as T = 10 2 k B T/e. Note that, 
T = 3.0 corresponds to room temperature. 

In the following sections, we study separately the effect of temperature and the 
role of intermolecular interactions by comparing systems with weak x-r bonds 
(e x . r = O.le) and broken x-x bonds (D = 0) to the original one (Tabled]). 



3 Static fingerprints of the glassy state 

We describe the static structure by calculating the distribution function of the 
polar angle 6 between the director n of the rod and the z axis (Fig. CP, as 

i M iVrod 

9( e ) = T7^Y,Y,tm J )-e), (5) 

MiV rod j=l i=1 

where M is the total number of time origins. A similar function g((p) can be 
constructed for the azimuthal angle <p in the xy-p\ane. 

The temperature dependence of the orientation distribution function g{9) is 
shown in Fig. HI At high temperatures (T = 1.9,3.0) a broad angular distri- 
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bution of rods orientation peaked at 9\ ~ 16° results from the sampling by 
the system large part of the available microstates. As temperature decreases, 
a second preferable orientation of the rods at 9 2 ~ 19° appears, indicating 
some structural rearrangement below T = 1.9 as seen at T = 0.9. At lower 
temperatures both peaks become more pronounced, until at T = 0.11 they 
merge into a "plateau" around 9\ and a peak at larger 9 appears. A crossover 
between high and low temperature structures seems to occur between 1.9 and 
0.9. This crossover is confirmed by looking at the azimuthal angle distribu- 
tion g(ip), shown in Fig. We see that while three peaks are present at all 
temperatures, the relatively sharp peak at <p = 225° for T < 0.9 determines a 
prefered orientation of the rods in the xy-p\a.iae. Consequently, both distribu- 
tion functions suggest the appearance of additional order at low temperatures 
below T = 1.9 as observed at T < 0.9. 

In order to single out the role of the xyloglucan network, we study two systems 
with broken x-x and with weak x-r bonds. The resulting distribution g{9) for 
the system with broken x-x bonds is rather similar to the original one (Fig.Hb). 
However, for the system with weak x-r bonds, g{9) has only one pronounced 
preferred orientation 9 2 (instead of being double-peaked), which suggests an 
easier reorientation of rods as for higher temperatures. The structure of the 
network of xyloglucan particles follows the structure of the rods as it is clearly 
illustrated in Fig. [2], where one can recognize the seven peaks in the distri- 
bution of distances within a triplet for all studied temperatures (for example, 
we show T = 3.0; 0.45). For low temperatures the distribution is sharper, im- 
plying that the network of xyloglucans is less flexible. From these results we 
can conclude that the network of xyloglucan changes the optimal structure of 
the system, inducing glassy features like a broad, or double peaked, angular 
distribution. 



4 Slow relaxation 

We study the dynamics of our system by calculating the mean-square dis- 
placement (MSD) and the orientational time correlation function of the rods. 
At all studied temperatures, we find that the slow relaxation processes are 
well described by power laws and logarithmic functions. In Fig. Owe show the 
time dependence of the MSD = (jr(t) — r(0)| 2 ^ of rod beads in log-log scale. 
This dependence can be described by a power law 

(|r(t)-r(0)| 2 ) = a-t 6 (6) 

with exponent b ~ 0.8 (instead of b = 2 characteristic for the ballistic regime) 
for the initial part (t < 20r) and with a temperature dependent b(T) for the 
slow regime (t > 30r) shown in the inset of Fig. The diffusive regime (6=1) 
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typical of liquids is never achieved in our simulations since, even at the highest 
temperature, the value of b is less than 0.5, signaling a slow glassy dynamics of 
the rods. Moreover, the appearance of a "plateau" between the two regimes for 
T = 0.9 corresponds to the emergence of another structurally arrested state, 
where the rods are trapped by their neighbours. This confirms that T = 0.9 
is below the crossover temperature between two glassy phases. By fitting b(T) 

(inset of Fig. [HD to a power-law we find b(T) oc \]~T. Using this approximation, 
one can estimate the temperature typical of diffusive dynamics as T ~ 12. We 
did not simulate such high temperature, where most xyloglucan bonds would 
be broken, since this regime is not interesting for cellulose. We found the 
following influence of xyloglucan network on the MSD at T = 0.45: weakening 
of x-r bonds changes only slightly the MSD, while broken x-x interactions 
change qualitatively the MSD behaviour and seem to suppress the subdiffusive 
dynamics. 

Another quantity that can be used to characterize the glassy behaviour is the 
relaxation of rod orientation given by the following correlation function of the 
directors n 

o M N Tod 

(i> 2 (n(t) ■ n(0))) = £ £ (mfe) • n^. + tj) - -, (7) 

ro j — 1 % — 1 

where P 2 stands for the Legendre polynomial. From Fig. [7] we notice that the 
above quantity has a behaviour similar to that of the MSD with two distinct 
relaxation processes: an initial fast exponential and slow logarithmic long-time 
decay. Thus, we fit these curves over the whole time interval with the following 
expression 

(P 2 (n(t) • n(0))) = A - c • log* +p ■ exp(-t/r). (8) 

The fitting parameter c that characterizes the slow dynamics is shown as a 
function of temperature in the inset of Fig. The curve c(T) saturates for 
T > 0.9, another sign of a structural rearrangement. Above this temperature, 
we find an increase of broken x-x bonds from 2% up to 6% at T = 3.0. The 
increased number of broken bonds at higher temperature makes the system 
less frustrated preventing further increase of c with temperature. Nevertheless, 
the slow logarithmic reorientation dynamics persists also above the crossover 
temperature. Together with the data on strongly subdiffusive translational 
dynamics at all temperatures, this suggest a crossover between two types of 
structurally arrested glassy states like the one observed experimentally in a 
copolymer mi cellar system [18j . 

The reorientation of rods depends also on the strength of the x-r interac- 
tions. Tenfold weakening of x-r interactions at T = 0.45 results in a change of 
P 2 (n(i) • n(0))) more than that due to an increase of temperature to T = 3.0 
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(Fig. [7]) but slow logarithmic dynamics still exists. Conversely, breaking of x-x 
bonds (D = 0) leads to a non decaying amorphous behaviour of the orientation 
of rods and eliminates the slow dynamics completely. 



5 Discussion and Conclusions 

We have studied by means of Molecular Dynamics simulations a model that 
captures the main structural features of cellulose. We have shown that, in a 
wide range of temperatures, a slow dynamics results from the competition be- 
tween the microfibrils interactions and the stretchable network of xyloglucan 
molecules. The slow dynamics is characterized by logarithmic time depen- 
dence of orientational correlators and by strongly subdiffusive dynamics of 
translational motion. These two features are robust and were observed for 
all temperatures investigated, around and below room temperature. More- 
over, we have shown that weakening of the microfibril-xyloglucan interactions 
preserves the slow dynamics but influences the time scale of diffusion and re- 
orentation. These findings are compatible with the observation of logarithmic 
creep motion in plant cell walls, also when in presence of LTP that weakens 
the hydrogen bonds in the network [3]. The loosening effect of LTP, rather 
than temperature, is thought to be one of the mechanisms that makes the 
extension of plant cell walls possible [5j. Indeed also in our model, variations 
of temperature around room temperature do not affect noticeably the slow 
dynamics in the way a weakening of x-r interaction does. 

The complexity of our model produces an additional structural transition at 
low temperature. We find that a crossover beetween two glassy states, char- 
acterized by different angular distributions and parameters of slow dynamics. 
A similar glass-glass transition was found in copolymer micellar system |18j . 
We do not pursue further the study of the nature of this transition, that in 
our model occurs at low temperatures, not relevant for biological processes. 
However, the possibility of more than one structurally arrested state can be 
relevant for other related biopolymer networks. 
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Table 1 



Parameters of interaction potentials 


Bond 


Potential 


Energy 


Length 


Cut-off radius i? CjS 


b-b 


Ub, ^LJ.b-b 


A = 2500, £b-b = 10 


(T b -b = 10 


-R c ,b-b = 15 


r-r 




e r _ r = 0.001 


<7 r _ r = 4 




x-r 




S X _ T = 1 = £ 


0"x-r = 1 = CT 


Rc,x-r — 2.3 


x-x triplets 


U M 


D = 0.05 


r eq G [20; 60] 


6 = 0.25 a' 1 (Eq.HJ) 


x-x all 


L^LJ.x-x 


e x _ x = 0.1 


0"x-x = 1 


-Rc,x-x =1.12 




Fig. 1. Schematic representation of our model: each rod are formed by seven beads 
(grey spheres). Each x-particle (black spheres) interacts with the rod beads and 
with other two x-particles (see text) as indicated by the dashed-dotted wavy lines. 
The ratio between rod beads and x-particles is 7:3. The rods are closely packed, i.e. 
the equilibrium length of rod L = 67.3a is significantly larger than the equilibrium 
space between them h = 4.5a. The rod director n is also indicated in cartesian and 
polar coordinates. 
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40 

length/ <7 

Fig. 2. Bond length distribution of connected x-x triplets. Dotted line: initial sam- 
ple providing the r eq -distribution (see text); dashed and solid lines: equilibrated 
system around room temperature (T = 3.0) and at low temperature (T = 0.45), 
respectively. 




Fig. 3. Snapshot of the system at T = 0.45 in two representations, (a) Yellow balls 
and sticks: rods; red balls: xyloglucan particles. One can distinguish two layers of 
rods roughly oriented along the z direction. For illustration purposes the x and y 
sides of the box are multiplyed by five, (b) The simulation box is shown with real 
relative dimensions. Only the x-particles (red balls) and a few of the bonds (black 
lines) within triplets are shown. 
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Fig. 4. Polar angle distribution g(Q) of rod orientation for the indicated reduced 
temperatures. The dotted vertical lines at 0\ ~ 16° and 82 ~ 19° indicate the two 
preferable orientations of the rods. In panel c (T = 0.45) the result for the system 
with weak x-r bonds (dashed line) and broken x-x bonds (dotted line) are shown 
as well. 




Fig. 5. Azimuthal angle distribution g(<p) of rod orientation for the indicated reduced 
temperatures. Notice the increase of the peak at ip = 225° for T < 0.9. 
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Fig. 6. MSD of rods in log-log scale. The initial part for all studied reduced temper- 
atures is described by a power law with exponent b ~ 0.8. The long-time behaviour 
of the MSD is described by a power law with temperature dependent exponent b(T) 
given in the inset. Notice, that for all temperatures the dynamics is strongly subdif- 
fusive (b < 1). The two time regimes are separated by a "plateau" at T = 0.9; 0.45. 
The MSD for a system with broken x-x bonds is not well described by a power 
law for t > 30r, while weakening of x-r bonds changes only slightly the MSD at 
f = 0.45. 




Fig. 7. Orientational time correlation function of the rod director n (Eq. [7|). The 
initial decay is exponential, whereas the long-time decay shows slow logarithmic 
relaxation. The breaking of x-x bonds (dotted line) at T = 0.45 eliminates the slow 
dynamics, while the weakening of x-r interactions (short-dashed line) leads to a 
higher value of c ~ 0.017 (Eq. [8]). Inset: the rate of change in rods orientation c(T). 
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